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Abstract 

We consider linear and nonlinear hyperbolic SPDEs with mixed derivatives with ad¬ 
ditive space-time Gaussian white noise of the form Yxt = F(Y) + crWxt- Such equations, 
which transform to linear and nonlinear wave equations, including Klein-Gordon, Liouville’s 
and the sine-Gordon equation, are related to what Zimmerman (1972) called a diffusion 
equation. An explicit numerical scheme is employed in both deterministic and stochastic 
examples. The scheme is checked for accuracy against known exact analytical solutions for 
deterministic equations. In the stochastic case with F — 0, solutions yield sample paths 
for the Brownian sheet whose statistics match well exact values. Generally the bound¬ 
ary conditions are chosen to be initial values Y{x, 0) and boundary values Y (0, t) on the 
quarter-plane or subsets thereof, which have been shown to lead to existence and uniqueness 
of solutions. For the linear case solutions are compared at various grid sizes and wave-like 
solutions were found, with and without noise, for non-zero initial and boundary conditions. 
Surprisingly, wave-like structures seemed to emerge with zero initial and bondary condi¬ 
tions and purely noise source terms with no signal. Equations considered with nonlinear 
F included quadratic and cubic together with the sine-Gordon equation. For the latter, 
wave-like structures were apparent with a < 0.25 but they tended to be shattered at larger 
values of a. Previous work on stochastic sine-Gordon equations is briehy reviewed. 
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1 Introduction 


There have been many recent articles involving the applications of stochastic partial differential 
equations. These include parabolic equations which have enjoyed widespread attention in 
neurobiological applications (Boulakia et ah, 2014; Dorsek et al, 2103; Faugeras and MacLaurin, 
2014; Khoshnevisan and Kim, 2015; Petterson et ah, 2014; Stannat, 2013; Tuckwell, 2013a, 
2013b) and to a lesser extent hyperbolic equations (Hajek, 1982). These applications have 
often employed two-parameter Wiener processes, or Brownian motion, {hF(x,t),x G X, t G T}, 
with mean zero and covariance Cov[lT(x, s), hF(y, t)] = min(rE, y) min(s, t) where X and T are 
sub-intervals of R or R'^ (or their formal derivatives, space-time white noise {w{x,t)}). Such 
two-parameter processes first appeared over 50 years ago in works by Kitagawa (1951), Cencov, 
(1956 ) and Yeh (1960). 

Zimmerman (1972) constructed a stochastic integral with respect to W and obtained a 
solution to what that author called a diffusion equation, 

AY{x,t) = Y(x,t) - Y(x,0) - Y(0,t)+Y(0,0) 

pt px pt px 

= / / m{u,v,Y{u,v))dudv + / / a{u,v,Y{u,v))dW{u,v) (1) 

Jo Jo Jo Jo 

where m and cr are Baire functions satisfying suitable growth conditions. It was shown, inter 
alia, that under the stated conditions, solutions Y had sample functions which were almost all 
continuous and were uniquely determined. 

Yeh (1981) considered solutions of a stochastic differential equation on H = [0,oo] x [0, oo] 
written as 

dY (x, t) = m{x, t, Y)dxdt + a{x, t, Y)dW (x, t) (2) 

with a boundary condition 

Y{x,t) = Z{x,t), ioi{x,t) £ dD (3) 

where dD is the boundary of D and Z is a random process with continuous sample functions and 
locally bounded second moment on dD. It was proved that under suitable growth conditions 
on m and a that a strong solution existed with pathwise uniqueness. Interestingly, it was stated 
that the system is non-Markov (Yeh, 1981). 

It seems that the stochastic equations (1) and (2) have not been employed as written in any 
biological, engineering or physical science modeling although we shall see that they transform 
to some well-known equations of mathematical physics. 


2 A simple stochastic PDE with mixed derivatives 

The integral equation (1) or the differential equation (2) can be formally recast as the stochastic 
partial differential equation 

d'^W 

— + — . (4) 

This is a hyperbolic equation which is similar to the general class of stochastic partial differential 
equations considered by Hajek (1982) 

d'^Y dY dY 
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where w is a space-time white noise with mean zero and covariance function 


Cov[rc(a:, s),w{y, t)] = 5{x — y)S{s — t). 


(6) 


In Hajek’s equation the drift and diffusion coefficients depend only on Y. Hajek also proved 
existence and uniqueness of solutions and described a Stratonovich type calculus. 

Identifying t as a time parameter and x as a spatial parameter, in the rest of this article 
we are concerned with temporally and spatially homogeneous equations with the structure 


d^Y 

dxdt 


F{Y)+a{Y) 


d'^W 

dxdt' 


( 7 ) 


which is obtained from Hajek’s equation by putting a = 0 and substituting F and a for h and 
c, respectively. 


3 Linear examples with space-time white noise 


A simple equation with the structure of Eq. (7) is the linear stochastic PDE with F(Y) = aY 
and (j(Y) = cj, where a and a are real constants. Thus 


d‘^Y _ d'^W 

dxdt ^ dxdt' 

or, using subscript notation for partial differentiation. 


( 8 ) 


Yxt = aY + aWxt- 


(9) 


Generally, in physical or natural sciences, it would be of interest to find solutions on [0, x/] x 
[0,tj] where Xf < oo and tf <oo with given values of Y, possibly random, on the boundaries. 

When fj = 0, y is deterministic and an exact general solution may be written as the linear 
combination 

y (x, t) = Cl exp(ax) exp(t) -|- C 2 exp(x) exp(Q;t) (10) 

Cl and C 2 are arbitrary constants. Note that an arbitrary constant cannot be added to this 
general solution and the boundary conditions giving values of Y (0, t) and Y (x, 0) are, if this 
form of the solution is used, being determined by the values of ci and C 2 , are not arbitrary. 


3.1 Numerical scheme with cr = 0 

To solve Eq. (8) without noise numerically on [0,x/] x discretize x and t with steps of 

Ax and At, so that Xi = {i — l)Ax, f = 1, 2,... , n^; -|- 1 and tj = (j — l)At,j = 1, 2,..., -|- 1, 

so that Ax = Xfjux and At = tf/nt- Then approximating Y{x,t) at the grid point {xi,tj) by 
Yij and using the definition of second order partial derivative leads to the scheme, 

Yij = Yi_ij + Yij_i + AxAtF(yi_ij_i) - Tj-ij-i. (11) 

Supposing the initial values (IC) at t = 0 are 

y(x,0) = /(x),0 < X < Xf < oo, (12) 

and the boundary values (BC) at x = 0 are 

y(0,t) = g{t),0 <t<tf< oo, (13) 
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with the necessary consistency condition /(O) = ^(O). Then the values of Tij on the boundaries 
are = f{xi),i = 1,..., re^; + 1 and Yij = g{tj),j = 1,..., re^ + 1. Such forms of initial and 
boundary values for linear and nonlinear PDEs of the form Y^t = F(Y) have been employed 
by Fokas (1997), Leon and Spire (2001), Leon (2003) and Pelloni (2005), mainly for the sine- 
Gordon equation. Pelloni (2005) pointed out that the use of such IC/BC gives a well-posed 
problem and no additional constraints are required or necessary to define a solution. In fact, 
the solution at (xo,to) depends only on solution values at x < xq and t < to so that one 
may continue to integrate the PDF indefinitely away from the IC/BC conditions (12) and 
(13). Given liq, Vi and Yi ^2 enables the calculation of Vi. Similarly, given TijVj and 12,1 
enables the calculation of Y 2 j,yj. Proceeding in this fashion, Yij is determined Vi, j from the 
initial/boundary values. 


C,=1, C2=0, a=1, a=0 c,=0.5, C2=3, a=-2, o=0 



Figure 1; The top two panels show the numerical solution of the linear PDF (8) without noise 
and with the parameters stated on the figure. The lower two panels show sample paths obtained 
by numerical solution of (8) with the same parameters but including noise with a = 0.5 on the 
left and a = 0.1 on the right. 

Fxact deterministic solutions, that is, for a = 0, were checked against the solutions of the 
form of Fqu. (10). For a first example, putting ci = 1, C 2 = 0 and a = 1 on [0,2] x [0,2] with 
nx = nt = 100 or Ax = At = 0.02 gave Y (2, 2) = 53.290 which compares with the exact value 
of = 54.598. However, with Ux = nt = 500 or Ax = At = 0.004, the numerical value was 
54.33, being much closer to the exact solution. The numerical solution with Ux = nt = 100 is 
shown as the surface in the top left part of Figure 1. As a second example, take ci = 0.5, C 2 = 3 
with a = —2 on [0,1] x [0,1] with nx = nt = 100 or Ax = At = 0.01. The numerical scheme 
gave T(l, 1) = 1.2728 compared with the exact value 3.5e“^ 1.2876. With nx = nt = 500 or 
Ax = At = 0.002, the numerical scheme gave Y (1,1) = 1.2846. The maximum value of Y over 
[0,1] X [0,1] was 8.2225 for both grid sizes. The numerical solution is shown for the grid with 
nx = nt = 100 in the top right part of Figure 1. 

As a further test of the numerical scheme, the same equation was considered on [0,2] x [0, 3] 
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with parameters a = 0.5, ci = 0.5, C 2 = 1. For a grid with Ux = nt = 500, the maximum 
relative error was 0.0035, and with Ux = 700, nt = 800 it was 0.0023. For a still finer grid 
Ux = 1000, nt = 1200 the maximum relative error was 0.0016. For this latter case the exact 
solution 

Y{x, t) = 0.5 exp(0.5x) exp(t) + exp(x) exp(0.5t) (14) 

evaluated at grid points and called Z{i,j) is shown in the top part of Figure 2, the bottom part 
giving the relative error against the exact solution. As the mesh becomes finer, the relative 
error diminishes at all grid points. The maximum value of 0.16% for the case shown points to 
the accuracy of the numerical scheme. 


3.2 Stochastic case 


To simulate the additional white noise term in Equation (9), we approximate as in the 

numerical solution of the spatial Fitzhugh-Nagumo equation (Tuckwell, 2008) with space-time 
white noise, where results were checked against analytical solutions for the moments. Thus at 




d‘^W _ 1 

dxdt 


(15) 


where where the Ajj are independent standard (zero mean, unit variance) normal random 
variables, which will be generated by computer random number generator. Thus the discretized 
version of (8) is 


Yi^j = Tj-ij + AxAtF(Fi_ij_i) - + ctV AxAtNi- 


A sample path with the inclusion of noise with cj = 2 for the first deterministic example 
(ci = 1, C 2 = 0, a = 1) is shown in the bottom left panel of Figure 1 (below the deterministic 
solution). Here again nx = nt = 100 and Xf = tf = 2. For 100 trials, the standard deviation 
of y(2,2) was 12.75 and the mean of y(2,2) was 54.66, compared to the no noise numerical 
value of 53.29 with the same grid size and the exact value of 54.60. 

In the bottom right panel of Figure 1 is shown a sample path, with cj = 2, for the second 
deterministic example (ci = 0.5, C 2 = 3,a = —2), again on [0,1] x [0,1] with nx = nt = 100. 
For 2 sets of 100 trials, the standard deviation of Y{1, 1) was 1.406 and 1.265, with means of 
0.976 and 1.20, respectively, the latter results being considerably less than the exact mean of 
1.2876. When the grid was made finer with nx = nt = 200, the mean was closer to the exact 
result at 1.26 and the standard deviation was 1.477 (more details below). 


3.3 Brownian sheet sample paths 


It was of interest to use the numerical scheme of Eq. (16) with F = 0 together with boundary 
values of 0 to generate sample paths for the two-parameter Wiener process or Brownian sheet 
(Weiner, 1975; Adler, 1978; Koshnevisan, 2001). This was done with u = 3 on [0,1] x [0,1] and 
on [0, 2] X [0, 2] for two different mesh sizes, nx = nt = 100, nx = nt = 400, with 500 and 200 
trials respectively. Examples of the sample paths are shown in Figure 2, the top two panels 
being for W on [0,1] x [0,1] and the bottom two on [0,2] x [0,2]. 

Table 1 gives the values of the mean E[IT {xf,tf)] and the standard deviation SD[W {xf,tf)] 
of W{xf,tf) from the simulations. The value of the exact mean is E[IT(x/,t/)] = 0 and that 
of the standard deviation is SD[W{xf,tf)] = xja, being either 3 or 6. The approximate 95% 
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confidence intervals for the mean are given in column 5. For both sets of trials on [0,1] x [0,1] 
the sample mean is well within the 95% confidence limits, but in the first run with 200 trials 
with the finer mesh on [0, 2] x [0, 2], the sample mean is just outside the 95% confidence interval. 
In a second set of 200 trials, the sample mean was well within the confidence limits. 


Table 1: Numerical results for simulation of Brownian sheet: in all cases <7 = 3 


rix = nt 

No. Trials 

II 

E[Wixf,tf)] 

95% C.I. for E 

SD[Wixf,2tf)] 

100 

500 

1 

-0.165 

±0.263 

3.104 

400 

200 

1 

0.237 

±0.416 

3.296 

100 

500 

2 

0.041 

±0.526 

5.803 

400 

200 

2 

-0.893 

±0.832 

5.867 

400 

200 

2 

0.310 

±0.832 

6.589 


n„=n,=100 n„=n,=400 



0 1 



n„=n,=400 



0 0 


Figure 2: Four simulated sample paths of a Brownian sheet W generated by applying Eq. (16) 
with F = 0 and a = 3. The top two panels are on [0,1] x [0,1] whereas the bottom two are 
on [0, 2] X [0, 2]. The grid is finer (400 by 400) in the left hand cases, being 100 by 100 in the 
right hand cases. Boundary values W = 0 are applied on the x and t axes. 


3.4 Further numerical results for the linear stochastic Equation (8) 

Considering further the SPDE (8) with noise, it is of interest to see how the results might 
depend on grid size and number of trials. To this end, with ci = 1,C2 = 0 and a = 1, several 
runs were done on [0, 2] x [0, 2] as described in Table 2 with the stated results for the mean and 
standard deviation of T(2,2), being the maximum point in the deterministic case. In general 
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it is apparent that the value of Y{2, 2) is underestimated in all the deterministic cases, though 
the numerical value seems to approach the exact value as the grid becomes finer. Thus the 
error with a grid of 100 x 100 gives an error of 2.39 % whereas a grid of 1000 x 1000 gives an 
error of 0.25 %. 

As it is apparent (and see also Section 3.5), by taking expectations and integrating directly 
in Equ. (8), that the mean of Y with the additive white noise is the same as the deterministic 
value, then the same remarks apply to the mean in the stochastic case, as the difference between 
the exact deterministic value and the mean with noise at a = 0.1 drops (monotonically for these 
simulations) from 2.36 % for a grid of 100 x 100 to 0.63 % for a grid of 400 x 400. From the 
results of the last 4 rows in Table 2 where a = 0.5, there is no consistent pattern of changes in 
the mean relative to the deterministic value as the number of trials increases, though the set 
of results is very small. 


Table 2; Results for Equation (8) with a = 1 for various grid sizes and noise levels 


Run 

a 

Ux = nt 

E[Y(2,2)] 

SD[Y(2,2)] 

Exact 

0 

- 

54.598 

0 

Numerical deterministic 

0 

100 

53.290 

0 


0 

200 

53.936 

0 


0 

500 

54.331 

0 


0 

1000 

54.464 

0 

Numerical SPDE 





No. trials 





200 

0.1 

100 

53.310 

0.601 

200 

0.1 

200 

53.910 

0.640 

200 

0.1 

300 

54.133 

0.657 

200 

0.1 

400 

54.253 

0.645 

100 

0.5 

200 

53.787 

3.112 

200 

0.5 

200 

53.520 

3.172 

300 

0.5 

200 

54.166 

3.129 

400 

0.5 

200 

53.918 

3.154 


Further results were obtained for the stochastic PDE (8), integrated with various noise 
levels a from 0 to 3, using the same parameters as above but that in every case Ux = nt = 200 
and 200 trials. When a = 2, for example, the mean and standard deviation of Y (2,:) were 
plotted against time for t up to 2. The agreement between the mean and the exact deterministic 
value Y{2,:) = exp(2)exp(t) was excellent. The standard deviation grew in an approximately 
linear fashion over the same time interval. When the standard deviation of the solution value 
Y (2, 2) is plotted against noise level, with all other parameters fixed, the growth of the standard 
deviation SD[Y (2, 2)] is about linear with a and in fact a good fit to the curve is given by 

SD[Y{2,2)]^6a, (17) 

which can be compared with the value 2cj for the Brownian sheet. 



3.5 Other boundary conditions: waves and purely noise generated waves 

In the remainder of this section we further consider mumerical integration of the linear equation 

= aY + 13 + aW^f (18) 

Because the noise term contributes zero to the mean we expect that Y = E\Y] will satisfy the 
same equation as the deterministic solution, so that 

= aY+ 13. (19) 

With the parameters a = — 1 and /3 = 0, the veracity of this claim was tested on the square 
[0, 5] X [0,5] with the initial/boundary conditions 

y(x, 0) = 1, 0 < X < 5, y(o, t) = 1, 0 < t < 5, ( 20 ) 

and with a grid of 500 x 500. In Figure 3 are shown results for the deterministic solution (red 
surface) and the means for 50 trials with a = 0.05 (blue surface) and with a = 0.1 (green 
surface). It can be seen that the three surfaces are almost identical so that the mean can, with 
suitable geometry, mesh sizes and parameter sets be found as the deterministic solution. This 
also provides a heuristic test of the accuracy of a numerical scheme. 


Y(x,t) 

and 

E[Y(x,t)] 



Red:^=0 
Blue:a=0.05 


Figure 3: For Y^t = —Y + aWxt the deterministic solution is shown in red and the means for 
the stochastic cases a = 0.05 and u = 0.1 are shown in blue and green respectively. The three 
surfaces are practically the same. Initial/boundary values both unity. 


3.5.1 Previous boundary condition 

With boundary conditions of the form of (10), wave-like solutions were not apparent with the 
parameters a = —1, ci = 0.5, C 2 = 1, cr = 0. This was the case no matter how large a space- 
time interval was considered, up to [0, 20] x [0, 20] with a grid 500 x 500. Furthermore, the 
numerical solution agreed precisely with the analytical solution. 
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Figure 4: Numerical solution of Eq. (18) for various grid sizes without noise, cr = 0, and with 
IC/BC y(a:,0) = 1,0 < x < 30 and Y{0,t) = 1,0 < t < 30. Parameters a = —1, /3 = 0. A, 
grid 300 x 300, B, 700 x 700, C, 1100 x 1100. 


3.5.2 Waves with non-zero boundary conditions 

With the same parameters as in the previous example but with the initial/boundary conditions 

Y{x, 0) = 1,0 < X < 30, y(0, t) = 1,0 < t < 30, (21) 

with no noise {a = 0), waves of a sinusoidal shape form as depicted in the three panels of Figure 

4. Solutions were obtained for meshes with = n* = 300, (left), Ux = nt = 700, (middle), 
and Ux = nt = 1100, (right). The amplitude of the waves diminishes as the grid becomes finer, 
presumably reflecting greater accuracy of the solution. 

This was examined further by computing the solutions at three different mesh sizes (500 x 
500, 1000 X 1000 and 2000 x 2000) on a larger area, [0,40] x [0,40]. Results are shown in Figure 

5, where y(x,40) is plotted against x. In each case there are 12 peaks but as the grid becomes 
finer, and presumably the solution more accurate, the amplitude of the waves gets smaller and 
seems to attain a constant value, along with a fairly constant spatial and temporal frequency. 

The appearance of wave-like solutions is not surprising since under the coordinate transfor¬ 
mation 

C = x + t (22) 

Tj = X — t (23) 

the equation 

Yxt = F{Y) (24) 

becomes the general wave equation 

% - Ym = -F{Y) (25) 

which, in the case of a nonlinear F is called a nonlinear Klein-Gordon equation or with linear 
F simply the Klein-Gordon equation of mathematical physics. The linear case is also a small 
amplitude approximation to the sine-Gordon equation considered below. 

When 2-parameter white noise with amplitude a = 0.5 was present, again with a = — 1, /3 = 
0, the mean and variance were obtained numerically for 50 trials on [0,30] x [0, 30]. The result 
for the mean E\Y{x,t)\ is plotted in Figure 6, where again the grid is finer from left to right (left 
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Figure 5: Y(x,tf) versus x with xj = 40 and tj- = 40 for the solution without noise shown in 
Figure 5, but with meshes of various sizes. Top picture, Ux = nt = 500; middle, Ux = nt = 1000, 
and bottom Hx = nt = 2000. 


300 X 300, middle 500 x 500, right 700 x 700). In the mean the waves are still discernible and the 
amplitude again decreases as the grid becomes finer. There are 9 wavefronts discernible in the 
right hand part of Figure 6 (stochastic case) and in the middle part of Figure 4 (deterministic 
case), for which the grids are both 700 x 700. When these two plots are rotated in order to 
make a comparison clearer, for the stochastic case, the wave amplitude of E\Y (x, t] increases 
as the wave progresses, but it is not known whether this effect is real or due to an accumulated 
numerical error for large x and t. 



Figure 6: Mean of 50 trials for the solution of Eq .(18) with y(x,0) = 1,0 < x < 30 and 
y(0, t) = 1,0 < t < 30 with noise of amplitude a = 0.5. Parameters a = —1, /3 = 0. Grids: A, 
nx = nt = 300, B, Ux = nt = 500, C, nx = nt = 700. 
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3.5.3 Wave-like behavior with noise and zero bonndary conditions 

Having seen that solutions of = aY+aWxt, with a = —1, a = 0.5 and with initial/boundary 
condition as in (20) or (21) gave wave-like solutions, it was decided to examine solutions with 
the zero boundary conditions, 

y(x,o) = o,< X < 2o,y(o,t) = o,o < t < 20. (20) 

With no noise the solution is not surprisingly identically zero for all x and t. 

With noise level set at u = 0.5 patterns emerge which suggest wavefront activity. The 
mean of solutions for 50 trials is shown in Figure 7 for a grid 900 x 900. The pattern seems to 
self-organize into quite large wave-like segments with a wavelength which is fairly constant. 

To the right of the plot of E[Y{x,t)] (Figure 7B) there is a plan view (looking in the 
direction of -z) of where 0 < a < max(£^[y(x, t)]) and lAiz) is the indicator 

function taking the value 1 if z G A and 0 otherwise. That is, only the top portion of the 
surface lying above z = a is shown. Such a plot emphasizes the wave-like nature of E\Y{x,t)\. 
The value of a = 0.15. (In the caption of Figure 7 is also given an upper cut-off value, but this 
is above the maximum which makes it effectively infinite.) To further highlight the wave-like 
nature of (the mean of) the solutions, in Figure 7C a sequence of approximate wavefronts is 
marked in the case of a grid of 500 x 500. The evidence for a wave-like structure seems quite 
convincing. 

The appearance of these putative wave-like structures is unexpected. There is no “signal” 
because /3 = 0 and the initial/boundary conditions for Y{x,t) are identically zero. Hence the 
source of the wave-like structures is purely space-time white noise acted upon by the dynamical 
system. This response can be compared to the Brownian sheet sample paths with the same 
initial/boundary conditions as depicted in Figure 2, where no distinct pattern of wavefronts is 
apparent. It seems that the operator L defined through casting the SPDE as 

IS 

induces wave-like solutions from scattered local changes in Y coming from the Brownian sheet. 
The phenomenon which has the appearance of producing organized clusters is reminiscent, with 
some imagination, of theories of the creation of the “universe from nothing” (Vilenkin, 1982, 
1984; Krauss, 2012). 
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Figure 7: Wave like structures with zero boundary conditions. A. Plot of expected value 
of Y{x,t) versus x and t for the solution of Eq. (18) with y(x, 0) = 0,0 < x < 20 and 
y(0, t) = 0, 0 < t < 20 with noise of amplitude a = 0.5. Parameters a = —1, /3 = 0. Grid: 
fix = nt = 900. B. On the Plan view with E\Y{x,t)] cut off below a certain positive value 
to emphasize the wave-like structures, 0.15 < E\Y] < 0.706 = max(£'[y]). C. Plan view 
for a coarser grid, 500 x 500, with inserted curves which have been visually estimated to fit 
approximately to hypothesized wavefronts 

4 Simple nonlinear equations 

The above SPDEs have all been linear. In this section we briefly consider two nonlinear source 
functions E which are common in biological and physical models. 

4.1 Quadratic F 

The classical example of a PDE with a quadratic source function is the parabolic equation 

Yt = Yxx + kY{l - Y) (28) 

which was studied in the setting of spatial patterns of gene frequency by Fisher (1937) and 
Kolmogorov et al. (1937). In the setting of the present article we consider the SPDE 

y,t = A:y(i - y) + aWxt (29) 

which transforms to a nonlinear Klein-Gordon equation. 

4.1.1 Deterministic example 

For the deterministic equation with cr = 0 numerical integration was performed with the 
initial/boundary conditions of the form of (26) but with Xf = tf = 30 and with the two 
values 0.1 and 0.25 rather than 0. The results consisted of waves, similar to those shown in 
Figure 6, which settled into an approximately sinusoidal form. Sections are illustrated in Figure 
8A where are shown plots of Y at the middle x-value versus t, that is, Y(t, 15), the blue and 
red curves being for the initial/boundary values of 0.1 and 0.25 respectively. The periods for 
both cases are similar but the amplitude is smaller for the larger initial/boundary condition. 
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Figure 8; A. Values as functions of time of the numerical solutions of Yxt = Y{1 — Y) on 
[0,30] X [0,30] at X = 15. The blue curve is for the initial/boundary value of 0.1, whereas 
the red curve is for the initial/boundary value of 0.25. B. A sample path for the SPDE (29) 
with a = 0.06 and IC/BC value 0.1. This may be compared with the mean shown in the next 
diagram. C and D. Mean (left) and standard deviation (right) obtained by numerical solution 
(50 trials) of Yxt = Y(1 — Y) + crWxt with a = 0.06 and initial/boundary value of 0.1. 


4.1.2 Stochastic example 

Numerical integration of the equation (29) was performed with many values of the parameters 
Xf = tf,nx = Tit and a, in conjunction with initial/boundary conditions 

V(x,0) = 0.1,0 < X < Xf,Y{0,t) = 0.1,0 <t<tf. (30) 

Well behaved solutions were only found for rather small values of a. 

The behavior of solutions on [0,10] x [0,10] with a grid Ux = nt = 500 was investigated for 
various values of a. 50 trials were performed for each parameter set. 

A representative sample path (surface) is shown in Figure 8B with a = 0.06. The cor¬ 
responding mean and standard deviation are given in Figures 8C and 8D, respectively. The 
sample path does not depart greatly from the mean but in others the second wave front has 
fractured (not shown). The wave-like behavior of the standard deviation roughly minics that 
of the mean, but there seem to be additional weak wavefronts. 

The most interesting aspect of solutions, however, is the extreme sensitivity to the magni¬ 
tude of the noise parameter a near a critical value. Table 3 gives maxima and minima for E\Y] 
over the region of integration starting at a = 0.05 and increasing to a = 0.06467150. At the 
smallest of these values, the mean remains bounded for the whole region with a minimum value 
of 0.1000 and a maximum of 1.3364. When a exceeds a critical value of (about) (Tc = 0.06467147 
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the minimum of E\Y] becomes infinitely negative, but when a = cTc the minimum of E\Y] is 
finite and positive at 0.0997. 


Table 3: Some solution properties for Yxt = Y{1 — Y) + aWxt 


a 

max E[Y(x,t)] 

min E[Y(x,t)] 

Path in Figure 15 

0.05 

1.3364 

0.1000 

blue dot-dash 

0.06 

1.3131 

0.0981 

blue solid 

0.063 

1.3078 

0.0992 

red solid 

0.064 

1.3017 

0.0995 

black solid 

0.0645 

1.3200 

0.0993 

green solid 

0.0646 

1.3025 

0.0989 

mauve solid 

0.06465 

1.3021 

0.0996 

cyan solid 

0.06467 

1,3196 

0.0996 

yellow solid 

0.064671 

1.3201 

0.1000 

blue dashed 

0.0646713 

1.3109 

0.0997 

red dashed 

0.0646714 

1.3258 

0.995 

black dashed 

0.06467145 

1.3215 

0.0995 

- 

0.06467147 

1.3140 

0.0997 

- 

0.06467148 

1.3318 

—oo 

- 

0.06467150 

1.3095 

—oo 

- 


The results for E\Y (5, t)] for various values of a from 0.05 to Uc are shown in Figure 15, 
with an expanded portion on the right. It is apparent that there is no systematic change in 
the behavior of the solutions as a increases towards the critical value. The 50 sample paths 
for the smallest value of a at which the mean became infinitely negative (namely 0.06467148) 
were analyzed. Divergence to large negative values occurred in just two of the sample paths. 



Figure 9: A. Plots of expectations E\Y(5,t)] from simulations of Yxt = Y(1 — Y) + crWxt on 
[0,10] X [0,10] for various a and initial/boundary value of 0.1 according to Table 3. B. An 
expanded version over a small subinterval. C. Ten sample paths at x = 4 for solutions of 
Yxt = Y{1 — Y) + aWxt on [0,10] x [0,10] for a = 0.1 and initial/boundary value of 0.1. Most 
sample paths stay positive in a narrow band but a few paths (in red), called singular paths, 
take on negative values and eventually tend to — oo. 
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Further analysis was made for the case a = 0.1 which is considerably greater than the critical 
value (Jc- Of ten trials, two resulted in large negative values of Y. Trials were stopped when 
the values of Y became less than -20 and these are shown in Figure 16, for x = 4 as a function 
of t (up to f ~ 4.5 ) along with the 8 trials which did not become large and negative. It is not 
known what causes some paths to persist in a downward trend, although once Y < 0 the source 
function Y[1 — Y) is negative. Apparently there were no cases in which paths exceeded unity. 
This property is reminiscent of the anomalous paths in stochastic Hodgkin-Huxley systems 
(Tuckwell, 2007, Figure 3). 

4.2 Cubic F 

The cubic source function is common in simplified modeling in Neuroscience and Cardiology 
where it is used to endow mathematical models with a resting state and a threshold. Usually 
the principal variable (representing a voltage) is coupled to a recovery variable which is not 
included here. Thus the following SPDE is considered briefly 

Y^t = kY{l-Y){Y-Yi)+aW,t (31) 

where k, Yi and a are constants. The values k = 4 and Yi = 0.5 will be used throughout. 

4.2.1 Deterministic example 

For the deterministic case with a = 0, the equation (31) was integrated on [0,30] x [0, 30] with 
a grid 2000 x 2000 and with initial {t = 0) and boundary (x = 0) values of 0.1 and 0.6. For 
both IC/BC values wave solutions formed and attained apparently constant amplitude and 
wavelength as x, t became larger. The values of solutions are plotted as functions of t at x = 15 
in Figure lOA. For both IC/BC values the plots reveal oscillations of about the same period 
but different amplitudes. For the smaller IC/BC value the maximum and minimum of Y{x,t) 
are 0.1000 and -0.0544, and oscillations are about the lower equilibrium point Y = 0. With 
the larger IC/BC value the maximum and minimum are 1.1569 and 0.6000 (the IC/BC) value. 
The mean value of Y is 0.98 as the oscillations are roughly centered on the upper critical point 
at y = 1. 

4.2.2 Stochastic example 

For the cubic source as above, the effects of the inclusion of noise with a = 0.05 with an IC/BC 
value of 0.6 on [0,15] x [0,15] with grid Ux = nt = 750 were investigated with 50 trials. For 
the 50 sample paths the minimum value of Y was 0.5996, just below the IC/BC value of 0.6, 
and the maximum value was 1.1411, substantially above the upper critical value of 1. 

For the smaller noise amplitude of a = 0.025, the mean and standard deviation of Y are 
plotted in Figures lOB and IOC and show that the wavefronts in E\Y (x, t)] and the standard 
deviation are quite clear with some diffuseness at large x and t. 

For the larger value of a = 0.05, a sample path is shown in Figure llA, revealing discernible 
but irregular wavefronts. Figures IIB and IIC give the mean and standard deviation of Y for 
this case. The standard deviation starts to climb dramatically at large x and t values, becoming 
so large that it flattens the wavefronts of E[y]. However, it is not known without much more 
detailed analysis whether this enlargement of the standard deviation is a real chance effect or 
due to a deficiency in the numerical integration scheme. 

In order to see if using an IC/BC value below the 2nd critical point of y = 0.5 could even¬ 
tually result in a noise induced transition to large amplitude waves, the SPDE was numerically 
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integrated with an IC/BC value of 0.45 on [0,15] x [0,15] with Ux = nt = 750 and a = 0.05. In 
50 trials the maximum value of Y was 0.5111 which indicated the absence of large amplitude 
waves, When the noise level was increased to a = 0.06, in 10 trials the maximum value of Y was 
0.4671, but examination of paths suggested that a large amplitude response might eventually 
occur. Hence, again with a = 0.06 the region of integration was increased to [0, 25] x [0, 25] 
and a grid with Ux = nt = 1000 was employed. In each of ten trials the maximum value of Y 
exceeded 1.4, well above the upper critical point, indicating that large amplitude responses had 
occurred. The maximum value of the mean was 0.9257 and its minimum value was -0.1435. 
Examination of the ten sample paths showed that after a few small amplitude wave-like seg¬ 
ments, a sudden transition seemed to occur to large amplitude wave-like structures. However, 
further analysis is required to ascertain whether these transitions are due to purely stochastic 
effects. 



Figure 10: A. Values as functions of time of the numerical solutions of Yxt = 4V(1 — V)(y — 0.5) 
on [0,30] X [0,30] at a: = 15. The blue curve is for the initial/boundary value of 0.1, whereas 
the red curve is for the initial/boundary value of 0.6. B and C. Mean and standard deviation 
of Y{x,t) with the noise parameter a = 0.025 and IC/BC value of 0.6. 



Figure 11: A. A sample path obtained by numerical integration of Yxt = YY (1 — Y){Y — 0.5) -|- 
crWxt on [0,15] x [0,15] with a = 0.05 and an IC/BC value of 0.6. B and C. The mean and 
standard deviation of Y{x,t) from 50 trials. 
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5 Sine-Gordon equation 

The sine-Gordon equation, mentioned above, has a rich and interesting history cutting across 
many scientific and mathematical disciplines. It evidently made its first appearance in Bour 
(1862) and Enneper (1870) in the context of pseudospherical surfaces and then in Frenkel and 
Kontorova (1939) in the theory of crystal dislocations. 

In the 1960s it gained attention as a model for elementary particles starting with the works 
of Skyrme (1961) and Perring and Skyrme (1962). Of interest were the soliton solutions called 
kink and anti-kink. By definition, solitons are a special type of solitary wave which maintain 
their form after a collision and thus are ascribed particle-like properties. The sine-Gordon 
equation in Euclidean coordinates, 


4*tt 4^ XX — siiicj), 


(32) 


has the well-known kink/anti-kink (single soliton) solutions (Marchesoni et ah, 1988; Guarcello 
et ah, 2015) 


(pK+,K {x,t) = Aaician 


exp 


±(x — Xo — ut) 

yjl — V? 


(33) 


where u is speed and denote kink and anti-kink. There are also oscillatory soliton solutions 
called breathers which can be stationary, with frequency tc < 1, 


(j)^^{x,t) = darctan 


\/l — sin cut 

uj cosh(x\/l — w^) 


(34) 


or moving. 

The sine-Gordon equation is also satisfied by quantities associated with magnetic fields 
at Josephson junctions between two superconductors (Josephson, 1965; Lebwohl and Stephen, 
1967; Scott and Johnson, 1969). More recently the sine-Gordon equation has been mentioned in 
biological settings such as DNA (Yomosa, 1984). Some historical review is contained in Braun 
and Frenkel (1998) with reference to the relation between discrete and continuous forms. 


5.1 Some previous studies of stochastic sine-Gordon equations 

One of the first studies of perturbations of sine-Gordon solutions was that of Joergensen et 
al. (1982) in the context of Josephson transmission lines. Included in the PDE was a term 
representing a bias current, a loss term and a thermal noise term. The latter was taken to be 
a 2-parameter Gaussian white noise of the type used throughout the present paper. A similar 
system was considered for sine-Gordon particles in Bergman et al. (1983) who put 

4>tt - (pxx = -sm(j) + I' - G(t)t (35) 

where J- represents external forces, including noise, and G is viscosity associated with the loss 
term. Marchesoni (1986) considered the same PDE but without the loss term in order to see 
how a random force field might influence kink motion. An assumption that the two-parameter 
white noise could be factored into the product of two single-parameter processes was made and 
a Langevin equation with a one-parameter noise was analyzed for the relativistic momentum 
of the soliton solution. This approach was extended in Marchesoni et al. (1988) to include the 
calculation of nucleation rates. 

Biller and Petruccione (1990) also considered a stochastic sine-Gordon equation of the same 
form (35) as Bergman et al. (1983), but with an additive noise term which was assumed to have 
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only time-dependence, being a zero-mean 1-parameter Gaussian white noise. They simulated 
the stochastic PDE and also performed a perturbation analsis in order to ascertain the effects 
of random perturbations on the center of mass position and speed of the soliton solution. 

That approach was not dissimilar to that in another early article (Pascual and Vazquez,1985) 
on the effects of random perturbations in the sine-Gordon equation. These authors numerically 
and analytically studied the inclusion of both additive and multiplicative white noise. Thus 
they put 

(pit - 4>xx =+ (36) 

where / is the perturbative term. Since they were only concerned with small noise they reduced 
the problem to seeing how random perturbations would affect the speed and position of the 
deterministic soliton solution (kink or antikink). However, the noise was single-parameter 
Gaussian white noise which was inserted into a set of ordinary differential equations. 

Technical analysis of stochastic sine-Gordon equations and its generalizations has been 
performed. Recently, such endeavours have been pursued energetically; for example by Hairer 
and Shen (2015) for 2 space dimensions, Anton et al. (2015) with two-parameter white noise on 
the unit square, Guarcello et al. (2015) with a numerical study of breathers in long Josephson 
junctions and Huang et al. (2015) with a study of stochastic bifurcation theory. 

5.2 Deterministic example 

The original form of the sine-Gordon equation is 

Y,t = smY (37) 

where x and t are sometimes called light-cone coordinates or null coordinates. This is in 
distinction to the Euclidean coordinates which have been more commonly used in physics 
and engineering applications where the more familiar wave equation is as in (25), or more 
customarily written as Yu — Y^x = — sinV. Note that some authors have considered — sinV 
on the right hand side of (37) (for example, Fokas, 1997; Leon and Spire, 2001; Leon, 2003; 
Pelloni, 2005). 

To numerically integrate (37) we use IG/BC conditions as in (12) and (13), and note that 
existence and uniqueness of solutions for this equation with IG/BC of this type on 0 < x < oo, 
0 < t < oo were proved by Angelova et al. (2008) as a Goursat problem. The starting point 
of their analysis was to write the solution as an integral equation obtained immediately from 
(37), 

f-x f-t 

Y{x,t) = f{x) + g{t) — c + / / sinY{u,v)dudv, (38) 

Jo Jo 

where c = /(O) = ^(0). As an example of a deterministic solution we integrate (37) on 
[0,15] X [0,15] with a grid 1000 x 1000 and with initial value f{x) = 1 and boundary value 
g{t) = 1. The computed solution is shown in Figure 12A. In addition we have integrated the 
PDE with — sinV on the right hand side using the same IG/BC conditions. This solution is 
shown in Figure 12B. A cross-sectional view of the two solutions at the half-way value of x = 15 
is shown in Figure 12C. Of interest in the expanded picture of Figure 12D is the smooth way 
solutions leave the boundary point at t = 0. 
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Figure 12: A. Numerical solution of the deterministic sine-Gordon equation Yxt = sinY on 
[0,15] X [0,15] with an IC/BC value of 1. B. The corresponding result for Yxt = — sinT. C. 
Numerically obtained solutions for Yxt = sinT versus time at the half-way point x = 15. D. 
Plot to t = 1 showing how Y leaves the boundary value of y(15,0) = /(15) = 1. 

5.3 Stochastic examples 

The sine-Gordon equation with additive two-parameter white noise in the form 

Yxt = sin F aWxt (39) 

was integrated by the numerical method outlined above in Eq. (16). In all the results to be 
discussed in this subsection the noise amplitude is cr = 0.1 and the initial/boundary conditions 
are f{x) = 1 and g{t) = 1. 

Figures 13A and 13B show the mean E\Y{x,t)\ and standard deviation SDEV\Y{x^t)] 
based on 50 trials on [0, 20] x [0, 20] with a grid 800 x 800. The mean does not seem to be very 
different from the solution shown in Figure 12A for the noise-free case, although the amplitude 
of the oscillations in the mean appears to diminish as the wave progresses. The standard 
deviation commences with quite small values at small x and t and begins to steadily grow in 
magnitude to reach values which are about 4 times larger at (20,20) than (2,2). However, the 
standard deviation maintains some of the wave-like character of the mean. 

In Figure 13C is shown a sample path for Y as in Equ. (32) again with a = O.I and the 
same IG/BC as before on [0,20] x [0,20]. The wave-like character of the underlying noise-free 
solution is discernible, which is not expected to be the case with larger values of a. This was 
borne out by increasing a, firstly to 0.25, and then to 0.5 - see Figure 14. 
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Figure 13: A, B. Mean and standard deviation for numerically integrated stochastic sine- 
Gordon equation Yxt = siny + aWxt with a = 0.1 on [0,20] x [0,20] with a grid 800 x 800. 
Initial/boundary conditions are unity. C. A sample path for the process as in A,B. 


5.3.1 Larger values of a 

A sample path for a = 0.1 on [0, 20] x [0, 20] with a grid 1600 x 1600 is shown in Figure 14A. 
For the larger value of a = 0.25 the sample path shown in Figure 14B has a coarse and ragged 
wave-like structure with several gaps. The values of Y are mainly positive and do not exceed 
about 6, which maximum is slightly larger than that for a = 0.1. Thus for a < 0.25 the original 
wave-like appearance, (on [0,20] x [0,20]) is roughly preserved. 



Figure 14: Sample paths for the stochastic sine-Gordon equation Yxt = sinT + (TWxt and a grid 
1600 X 1600 on [0, 20] x [0, 20] with cr = 0.1 (A), a = 0.25 (B) and a = 0.5 (C), Initial/boundary 
conditions are unity. 

With the even larger value of u = 0.5 the sample path shown in Figure 14C has a startlingly 
different character, but having a generally wave-like appearance. Y is mainly above zero 
commencing with a furrowed L-shaped region (blue) of relatively small values which gives 
way to larger and more haphazard fluctuations as {x,t) —)• (20,20). In another example (not 
shown), Y was mainly positive but less than 10, and with a noticeable deep trench which 
extended down to near -10. In the positive regions there were wave-like structures but they 
were not patterned on the deterministic solution nor the mean of the small noise case. Thus 
the wave-like structure was shattered by the imposed fluctuations. 
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5.3.2 Increase in number of grid points 

In order to see if making the grid finer has a significant effect on the numerically calculated 
solution properties, 50 trials were executed for the small noise case a = 0.1 of the sine-Gordon 
equation (39) with a grid of 400 x 400 and one of 800 x 800, both on [0,10] x [0,10]. The 
standard deviations for these two cases were not dissimilar, but values with the finer grid were 
somewhat smaller. Wave-like structure was apparent in both results but it seemed to be better 
defined with the finer grid. 

6 Discussion 

As it stands, the form of general equation considered throughout this article without noise term 
is 

Ya:t = F{Y). (40) 

The best-known example is the sine-Gordon equation with F = sinT which made its first 
appearance over 150 years ago. However, all nonlinear (and linear) wave equations of the form 

</>« - (pxx = G{4>) (41) 

can be transformed to the form of (40) by rotating (and rescaling) the coordinates as for 
example in (22) and (23). 

All the numerical solutions in this article were obtained with MATLAB on a PC which 
puts limits on the mesh/domain sizes through limited memory. The question of accuracy of 
the simulations for determining properties such as mean and variance for the stochastic cases 
involves the fine-ness of the grid, through the magnitudes of Ax and At, and the number of 
trials. The magnitude of the fluctuations relative to the solution for cj = 0 is probably also an 
important factor. The convergence properties of the simple Euler scheme in the present context 
are not known, but in the linear case, when the noise is additive one may compare the mean 
with the deterministic solution as elaborated on at the start of subsection 3.5. It is likely that 
for small values of a in the nonlinear examples, the criterion of approximate agreement of the 
mean and the deterministic solution can also be used to gauge the accuracy of the numerical 
scheme if the number of trials is large enough. Specifically, Figures 6 and 7 for the linear SPDE 
point to the grid’s not being fine enough, as the domain is large, being [0,30] x [0,30] in Figure 
6, and [0, 20] x [0, 20] in Figure 7. In Figure 6 the means are compared for various grid sizes 
but even with the finest grid of 700 x 700 and 50 trials the mean and deterministic solutions 
differ considerably. Further computations are needed to see whether the disparity diminishes 
at larger numbers of trials and greater fine-ness of grid. 

For the quadratic SPDE, as seen in Figure 8, agreement of Y and Ya^=o was close, but for 
the cubic SPDE, Figures 10 and 11 show good agreement only for the first few wavefronts. For 
the sine-Gordon equation agreement between Y and Yu=q is good, as seen in Figure 13 where 
the domain is [0, 20] x [0, 20] and the grid 800 x 800 giving Ax = At = 0.025 which apparently 
gives good convergence for 50 trials and a = 0.1. 

The boundary conditions employed in the numerical calculations have mainly been constant 
Y along X = 0 and t = 0, with a consistency requirement at (0,0). The wave equation (41) 
is usually endowed with initial conditions for u and for ut along with boundary conditions if 
on a finite space interval. The relation between solutions of (40) and (41) will be explored in 
a future article. Of much interest will be the determination, in particular for the sine-Gordon 
case, the nature of solutions in (40) which correspond to soliton solutions in (41). 
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Another classical PDE of the form of Equ. (40) is Liouville’s equation which can be written 
as (Kamran, 2002, who calls Yxt = F{Y) an E-Gordon equation, the case F = 0 being the wave 
equation) 

Yxt = (42) 

with associated wave equation 

Yu - Yxx = -e^. (43) 

It also has its origins in the analysis of the Gaussian curvature of a metric in differential 
geometry. (Note that there are two other completely different equations called Liouville’s; one 
in statistical mechanics and the other in quantum mechanics). Liouville (1853) gave a general 
explicit solution of equation (42), but it is difficult to consider numerical solutions of either 
the deterministic equation or its stochastic forms because suitable boundary conditions are not 
available. Preliminary investigations resulted in solutions which rapidly became unbounded. 

Finally, we note that equations of the form of (40) or more generally 

Yxt = F{Y,x,t) (44) 

can constitute a general form of growth model where the space-time rate of change of a biolog¬ 
ical, chemical or physical quantity depends on population size as a function of space and time. 
However, most of the applications of (40) can be, or have been, obtained by the transformation 
of a wave equation. 
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